# Piece of code used to convert the normalized spectrum to a flux calibrated
# spectrum. Written by Sergio Scarano Jr (scarano@astro.iag.usp.br) jul 2009.

import numpy
import scipy

# Taking in to account the geometry of the source

if extended == 'y':
    magint=mag+2.5*numpy.log10(pixsize**2)
else:
    magint=mag

# Getting the transmission curve for the chosen filter

if obsfilter == 'g':
    namefilt = os.path.join(_path,'filter_g.dat')
elif obsfilter == 'r':
    namefilt = os.path.join(_path,'filter_r.dat')
elif obsfilter == 'U':
    namefilt = os.path.join(_path,'filterU.dat')
elif obsfilter == 'B':
    namefilt = os.path.join(_path,'filterB.dat')   
elif obsfilter == 'V':
    namefilt = os.path.join(_path,'filterV.dat')   
elif obsfilter == 'R':
    namefilt = os.path.join(_path,'filterR.dat')
elif obsfilter == 'I':
    namefilt = os.path.join(_path,'filterI.dat')

filtertransm=open(namefilt)

wavefilt=[]
transmis=[]

for ss in filtertransm.readlines():
    wavefilt.append(float(ss.split()[0]))
    transmis.append(float(ss.split()[1].strip()))

wavefilt=scipy.array(wavefilt)
transmis=scipy.array(transmis)

# Reading the source spectrum (the wavelength is measured in Angmstrons and the flux
# is relative to that measured in 5500 A). They were originally taken from the GEMINI
# webpage and their references are:

   # Elliptical galaxy - a template spectrum (class T=-5, -4) covering the wavelength
   # range 22 - 10000 nm and taken from the Pegase Atlas of Galaxies by Fioc &
   # Rocca-Volmerange (1997). The spectrum has a sampling of 1nm at wavelengths
   # shortwards of 878nm and 20nm longwards.

   # Spiral (Sc) galaxy - a template spectrum (type Sc, class T=5) covering the
   # wavelength range 22 - 10000 nm and taken from the Pegase Atlas of Galaxies
   # by Fioc & Rocca-Volmerange (1997). The spectrum has a sampling of 1nm at
   # wavelengths shortwards of 878nm and 20nm longwards.

   # QSO - covers the (rest) wavelength range 80 - 855 nm and is taken from Vanden
   # Berk et al. (2001, AJ, 122, 549). The spectrum has spectral resolution ~1800
   # and a sampling of 0.1nm.

   # HII region (Orion) - covers the wavelength range 100 - 1100 nm and is taken
   # from the HST STIS exposure time calculator. The spectrum has a sampling of 0.05nm.

   # Planetary nebula (NGC7009) - covers the wavelength range 100 - 1100 nm and is
   # taken from the HST STIS exposure time calculator. The spectrum has a sampling
   # of 0.05nm.

def readSpec(file):
	return numpy.loadtxt(file,unpack=True)[:2]

def bbSpec(wave,Temp):
	h_cgs = 6.6260755e-27	
	c_cgs = 2.99792458e10
	k_cgs = 1.380658e-16
	wave_cgs = wave * 1.0e-8 # considero que wave esta em angstrom e passo para cm
	a = (2.0*h_cgs*c_cgs**2.0) / wave_cgs**5.0 
	b = 1.0 / ( np.exp( h_cgs*c_cgs / (wave_cgs*k_cgs*Temp) ) - 1.0 ) 

	return a*b
	
wave=[]#scipy.array(wave)
flux=[]#scipy.array(flux)


if spiral == 'x':
   spectrum=os.path.join(_path,'scGalaxy.dat')
   wave,flux = readSpec(spectrum)
if elliptical == 'x':
   spectrum=os.path.join(_path,'eGalaxy.dat')
   wave,flux = readSpec(spectrum)
if hiireg == 'x':
   spectrum=os.path.join(_path,'orion.dat')
   wave,flux = readSpec(spectrum)
if pn == 'x':
   spectrum=os.path.join(_path,'pn.dat')
   wave,flux = readSpec(spectrum)
if qso == 'x':
   spectrum=os.path.join(_path,'qso.dat')
   wave,flux = readSpec(spectrum)
if star == 'x':
   spectrum=os.path.join(_path,'star.dat')
   wave,flux = readSpec(spectrum)
if blackbody == 'x':
   spectrum=os.path.join(_path,'star.dat') # Get wavelenght from stellar spectrum
   wave,flux = readSpec(spectrum)
   flux = bbSpec(wave,user_Tbb)

#import pylab as py
#py.plot(wave,flux)
#py.show()
#exit(0)
# Using a spline to represent the transmission curve inside the
# limits of the observed spectrum
from scipy.interpolate import splrep, splev
tck = splrep(wavefilt, transmis)
curve = splev(wave, tck)

ind=numpy.where((wave < wavefilt.min()) | (wave > wavefilt.max()))
curve[ind]=0.0


# Applying the transmission curve to the normalized spectrum

fluxfilt=flux*curve

# Selecting the flux at the magnitude m=0 [erg sec^-1 cm^-2 A-1]

if obsfilter == 'g':
    F0filt = 4.14E-09
    extcoef=0.14
elif obsfilter == 'r':
    F0filt = 3.00E-09
    extcoef=0.11
elif obsfilter == 'U':
    F0filt = 4.19E-09
    extcoef=0.60
elif obsfilter == 'B':
    F0filt = 6.60E-09  
    extcoef=0.40
elif obsfilter == 'V':
    F0filt = 3.61E-09
    extcoef=0.20
elif obsfilter == 'R':
    F0filt = 2.26E-09
    extcoef=0.10
elif obsfilter == 'I':
    F0filt = 1.23E-09
    extcoef=0.08

# Performing the calibration

execfile(os.path.join(_path,"slitseeing.py"))

#FluxSource=F0filt*10**(-0.4*magint) # [erg sec^-1 cm^-2 A-1]
FluxSource=(1-cloucov)*fracslit*F0filt*10**(-0.4*(magint+extcoef*acceptableairmass)) # [erg sec^-1 cm^-2 A-1]
FluxNorm=numpy.trapz(fluxfilt,wave)

factnorm=FluxSource/FluxNorm
calibflux=factnorm*flux

# Saving the calibrated spectrum into a file

table=numpy.zeros((wave.size,2))
table[:,0]=wave
table[:,1]=calibflux
filename=objname+'_Calib.dat'
#save(filename, table, fmt='%.6e', delimiter='    ')

#pylab.plot(wave,calibflux)
#show()

